#ifndef UTIL_H_
#define UTIL_H_
#include "./Constants.h"
//工具箱

namespace Util{
    double max_intensity(const double* Et_array, int array_size);
    double peak_electricField(const double* Et_array, int array_size);

    int Get_Nradius(double z0, double waist, double focus, double wavelength, double h_min, double h_max);
    // International System of Units    
    double Get_phase(complex data); 

    double Get_Pin(double waist0, double I0);
    double Get_Pcr(double n0, double k0, double n2);
    // double Get_zc(double zr, double Pin, double Pcr);
    double Get_zc(double zr, double focus, double Pin, double Pcr);
    //矩阵转行向量 rows,cols分别是矩阵array的行和列
    template <typename T>
    void vectorization(T** array, T* vector, int rows, int cols, int start_rows){
        for (int i = 0; i < rows; i++){
            int i_index = i + start_rows;
            for(int j = 0; j < cols; j++){
                int j_index = i * cols + j;
                vector[j_index] = array[i_index][j];
            }
        }
    }
    //行向量转矩阵+转置 rows,cols分别是矩阵array的行和列
    template <typename T>
    void trans_reshape(T** array, T* vector, int rows, int cols, int start_cols){
        for (int i = 0; i < cols; i++){
            int i_index = i + start_cols;
            for(int j = 0; j < rows; j++){
                int j_index = i * rows + j;
                array[j][i_index] = vector[j_index];
            }
        }
    }

    //行向量转矩阵
    template <typename T>
    void reshape(T** array, T* vector, int rows, int cols, int start_rows){
        for (int i = 0; i < rows; i++){
            int i_index = i + start_rows;
            for(int j = 0; j < cols; j++){
                int j_index = i * cols + j;
                array[i_index][j] = vector[j_index];
            }
        }
    }
    //拆分:将Array_all中从start_index开始的size个元素拆分到array中
    template <typename T>
    void splid(T *array_all, T* array, int start_index, int size){
        for(int i = 0; i < size; i++){
            int i_index = i + start_index;
            array[i] = array_all[i_index]; 
        }
    }
    //合并:将array中的元素合并到array_all中从start_index开始的size个元素中
    template <typename T>
    void combine(T *array_all, T* array, int start_index, int size){
        for(int i = 0; i < size; i++){
            int i_index = i + start_index;
            array_all[i_index] = array[i]; 
        }
    }
    
}
double Util::max_intensity(const double* Et_array, int array_size){
    double maxI = 0;
    for (int i=0; i<array_size;i++)
    {
        double I = 0.5 * Constants::epsilon_0 * Constants::c * std::pow(Et_array[i],2);
        if (I > maxI)
        {
            maxI = I;
        }
    }
    return maxI;
}

double Util::peak_electricField(const double* Et_array, int array_size){
    double maxE = 0;
    for (int i=0; i<array_size;i++)
    {
        double E = fabs(Et_array[i]);
        if (E > maxE)
        {
            maxE = E;
        }
    }
    return maxE;
}
int Util::Get_Nradius(double z0, double waist, double focus, double wavelength,
                           double h_min, double h_max){ // e指数步长函数 获得nr
    double E0 = 1.;
    double E1 = 1e-6;
    double zr = PI * waist * waist / wavelength;
    double wz = waist * std::sqrt(1. + std::pow((z0 - focus) / zr, 2));
    double r_max = wz * std::sqrt(std::log(E0 * waist / (E1 * wz)));
    
    double a = 1. + (h_max - h_min) / (r_max);
    double b = h_max / h_min;
    double N = std::log(b) / std::log(a);
    // cout<<"r_amx = "<< r_max<<endl;
    // cout<< "Nr = " << N <<endl;
    return round(N);
}

double Util::Get_phase(complex data){
    double real = data.real();
    double imag = data.imag();
    double phase = std::atan2(imag,real);
    return phase;
}

double Util::Get_Pin(double waist0, double I0){
    double Pin = PI * std::pow(waist0,2) * I0 * 0.5;
    return Pin;
}
double Util::Get_Pcr(double n0, double k0, double n2){
    double Pcr = 3.77 * PI * n0/(2 * std::pow(k0,2) * n2);
    return Pcr;
}
// double Util::Get_zc(double zr, double Pin, double Pcr){
//     double bM = 0.367;
//     double aM = 0.852;
//     // aM = 0.892;
//     // double bM = 0.405;
//     double zc = (bM * zr)/std::sqrt(std::pow(std::sqrt(Pin/Pcr)-aM,2) - std::pow((1-aM),2));
//     return zc;
// }
double Util::Get_zc(double zr, double focus, double Pin, double Pcr){
    double bM = 1;
    double aM = 0.5;
    double theta = focus / zr;
    
    // aM = 0.892;
    // double bM = 0.405;
    // double zc =(zr + focus*focus/zr) / (std::sqrt((Pin/Pcr)-1) + theta);
    double zc =aM*(zr + focus*focus/zr)/ (std::sqrt((Pin/Pcr)-1) + bM*theta);
    // double zc =(zr + focus*focus/zr) / std::sqrt(std::pow(std::sqrt(Pin/Pcr)-1,2)) + theta;
    return zc;
}



#endif